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Abstract. We show analytically and numerically that the appearance of lumps 

and gaps in the distribution of n competing species along a niche axis is a robust 

phenomenon whenever the finiteness of the niche space is taken into account. 

In this case depending if the niche width of the species a is above or below a 

threshold Uc, which for large n coincides with — , there are two different regimes. 

, For a > ac the lumpy pattern emerges directly from the dominant eigenvector of 

■^^ , the competition matrix because its corresponding eigenvalue becomes negative. 

• ■ For tr < (Tc the lumpy pattern disappears. Furthermore, this clumping transition 

S ' exhibits critical slowing down as a is approached from above. We also find that 

the number of lumps of species vs. cr displays a stair-step structure. The positions 
of these steps are distributed according to a power-law. It is thus straightforward 
to predict the number of groups that can be packed along a niche axis and it 
coincides with field measurements for a wide range of the model parameters. 
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1. Introduction 

An important problem in ecology is how closely can species be packed in a natural 
environment [1]. A usual way to approach this issue is by considering the species 
distributed along a hypothetical one-dimensional niche axis [T]. To fix ideas one may 
consider the niche axis as a gradient that is related to the size of organisms. Each 
species i is represented by a normal distribution Pi{0— ^^P hi^ ^ /^j)'^/2f^ ] centered 
at i^i, corresponding to its average position ^ on this niche axis, and with a standard 
deviation a, which measures the width of its niche. The competition for finite resources 
among the n species can be described by a Lotka-Volterra competition model (LVCM): 

dN, N " 

-^ = ^»;^(^*-2^«*i^j)> (1) 

* i=i 

where Ni is the density of species i, r^ is its maximum per-capita growth rate, Ki is the 
carrying capacity of species i and the coefficients a^- is the competition coefficient of 
species j on species i. It seems natural to assume that the intensity of the interaction 
between two species i and j depends on how close they are along the niche axis. A 
measure of this is provided by the niche overlap, i.e. the overlapping between Pi(^) 
and Pj{C)- The competition coefficients a^ can be computed by the MacArthur and 
Levins overlap (MLO) formula |2j: 






(2) 



Recently Scheffer and van Nes [3] found by simulations that the combination of 
LVCM ID) plus MLO © yields long transients of lumpy distributions of species along 
the niche axis [For asymptotic times, the lumps are thinned out to single species unless 
a stabilizing mechanism/term is included, as it was shown in [3].] 

This phenomenon of spontaneous emergence of self-organized clusters of look-a- 
likes separated by gaps with no survivors was dubbed by the authors as self- organized 
similarity (SOS). It was recognized as an important new finding in an established 
model in ecology [4j [5] In addition, there is empirical evidence for self-organized 
coexistence of similar species in communities ranging from mammal [6] and bird 
communities [7] to lake plankton [8]. 

However, there has been some controversy on whether this lumpy distribution of 
species is indeed a robust result or rather depends strongly on details of the model, 
like the competition kernel [HI HH] ■ 

Here we show that the lumpy pattern is a robust phenomenon provided one 
takes into account the finiteness of the niche axis. Thus, truncation besides being 
a crucial assumption which guarantees clustering, allows the analytical computation 
of the eigenvalues and eigenvectors of the competition matrix A with elements Qij 
given by ^. Furthermore, we show that ultimately solving the linear problem is 
enough to get both the transient pattern -lumps and gaps between them- as well as 
the asymptotic equilibrium. The plan of this work is as follows: 

Since an analytic solution for realistic conditions - species randomly distributed 
along a finite and non periodic niche axis, each with a different per capita growth rate 
Ti and carrying capacity Ki - is not possible, we will consider in section 2 a series of 
simplifications. We get an analytic expression for the state of this simpler system, in 
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terms of the dominant eigenvector of A. It provides a qualitatively good description of 
the system for not too short times and becomes very good for asymptotic times. Part 
of the material of this section was presented in a previous short paper [TT] , but there 
are some important differences like considering a less rough approximation together 
with some steps better explained. 

In section 3 we show, using simulations, that all these simplifications do not 
destroy SOS: lumps and gaps remain in the case of a finite linear niche axis no matter 
if the niche is non periodic {i.e. it has borders), or the species are randomly distributed, 
or r and K changes from species to species. Indeed we go further and show that SOS 
occurs in niches of more than one dimensions or when interaction kernels different 
from the Gaussian kernel are considered. 

In section 4 we show that the prediction of the number of lumps as a function of 
a is in good agreement with measures in several ecosystems [T] , provided a is greater 
than a threshold value Uc. For this critical value it occurs a bifurcation which is 
responsible for the clumping transition. 

Section 5 is devoted to conclusions and to put our results in its proper perspective, 
addressing some general concerns about SOS and comparing with other different 
approaches. 

2. AN ANALYTICAL PROOF OF SELF-ORGANIZED SIMILARITY 
IN A SIMPLIFIED CASE 

We start by considering the following simplifications: 

51 - The n species are evenly distributed along a finite niche axis of length L ^ 1, i.e. 
Hi = {i-l)/n {i=l,...,n). 

52 - To avoid border effects, the niche is defined circular, i.e. periodic boundary 
conditions (PBC) are imposed. This is done by just taking the smallest of |/ii-/^j| and 
l-|/Xi-/ij| as the distance between the niche centers. 

53 - All species have the same per capita growth rate which we take equal to 1, r^ = 
1 for all i. 

54 - The carrying capacity K is also homogeneous, Ki = K for all i. 

Under the simplifying conditions S3 and S4 the system of equations ([T]) reduces 
to 

-^=Xi{l-^a.,jXj), (3) 

where Xi is the density of species i, normalized by its carrying capacity Ki (xi = 
N^/K,). 

An equilibrium of the system Qis specified by a set of densities x*, one for each 
species i, verifying: 

n 

A standard procedure to check the stability of this equilibrium is linear stability 
analysis. That is, to consider, initially small, disturbances yi{0) from the equilibrium 
values X* and study their fate yi{t) as the time grows. Let's take x* = x*\fi which, by 
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virtue of conditions SI and S2, is an exact equilibrium |^. The evolution equation for 
yi{t) can be written as 



dy, 
dt 



^ = -{x* + y.,{t)) Y^ a,jyj{t). 



(5) 



j=i 



Since the coefficients of the matrix A given by ^ are symmetric, in the eigenvector 
basis {vi}, it becomes diagonal with all its eigenvalues Xi real. Hence integrating 
equation ([5]) , and using that yi{0) is small, yi{t) can be approximated by 



y,(i) ~ y,(0)e-- 



'Xit 



(6) 



Thus, for asymptotic times, y becomes proportional to the dominant eigenvector v™, 
the one associated with the minimum eigenvalue of A, Am, *-e. 



y(i) oc e 



-X* Xmt^JJl 



v"^ (for large times). 



(7) 



We will show that, for a wide range of the parameters n and a, Xm{n^a) is in general 
negative (see below). Hence, from ([7]), y is amplified over time instead of decaying 
to zero (as it would happen in the case of a positive A„i). Therefore, for large times, 
from ([5]) we can express the time derivative of x as 



d-K 

Tt - --(^)^"-" 

and by integration we get the approximated solution given by 

X K. e~'^"^ * (for large times). 



(8) 



(9) 



Analytic expressions for the eigenvalues and eigenvectors of A are not known for the 
general case of random distributions of species on a niche axis with arbitrary boundary 
conditions. However, for the simpler case when the n species are evenly spaced along 
the niche axis, /ij — {j — l)/n (with the index j=l,...,n), and PBC (the simplifying 
conditions SI and 5"^) A becomes a matrix whose rows are cyclic permutations of the 
first one: 



Cl 



C2 
Cl 



C2 C3 



Cn — 1 ^n 

Cn-2 Cn-1 



Cl 



with Cj{n,a) — e^^^s^^ , where the tilde stands for (mod ^^^) implementing then 
PBC. For this case, the eigenvalues Xk and the components of the eigenvectors v'' [k 
= l,...,n) are given by pT^ : 



n 



Afe = ^c,(n,a)e*2^('=-i)^^ 






§ We later checked by simulations that all the derivations below are independent from the initial 
condition: the same results are obtained when starting from a completely random assignation of 
densities. 
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(10) 



and 



vf = n"^[cos(27r(fc-l)/^j) + sin(27r(fc-l)/Xj)] 



cos (MUzMz}! 

\ n 



27r(fc-l)(j-l) 



(11) 



Since the matrix A is symmetric Cj — Cn+2-j- Therefore, from (jlOp one can see that 
the eigenvalues occur in pairs: A^ ~ /^n-k, with the exception of Ai (and of A„/2+i if 
n is even). Furthermore, these paired eigenvalues can be expressed as 



n/2 

Afc = 2^Cj(n,cr)cos[27r(fc- l){i - l)/r 
i=2 



(12) 



Equation ()10p can be used to determine the index k = m that gives the minimal 
eigenvalue, for n and a given, \m{n,a) (as we have just seen, the index k = n-m+2 
produces the same value). The surface depicted in Fig. [T] corresponds to X„i{n,a) 
computed for a grid 2< n < 200, and 0.05 < a < 0.5. Notice that A„i is negative 
except for small values of a and becomes positive when n < 8. The substitution of 
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Figure 1. The minimal eigenvalue of A, Xm , determined from equation 1101 as 
a function of n and O". The black spot denotes the point n = 200 and a = 0.15. 
Inset: a zoom of Am vs. a for n = 200 in the interval 0.05 < cr < 0.1. 
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the dominant eigenvector v™, which from dill) has m-1 peaks and m-l valleys, into 
^ allows to predict the distribution of species for long enough times. 

The results we got were checked by numerical simulations. In these simulations 
the initial values for the Xi are random numbers between and 1. Then the system of 
differential equations (ODE) is integrated for a given final time. In Fig. [2]we compare 
this analytical approximation with simulations. For instance, if n = 200 and a — 0.15 
we get TO = 5 ( and m = 200-5+2 = 197), Am = 0.3938 and the components of v™ are 

given by W ^ sin[87r/ij] + i M cos[87r/ij]. Panel (A) of Figl2]if for t =1000 generations. 
The agreement is quite good and the quality of the agreement improves with time, 
until it becomes very good when the lumps are thinned to single lines as it is shown 
in panel (B j]J. This happens because we are not considering any lump stabilizing term 
like the one considered in [3]. Notice that ultimately the lumps and gaps coincide, 
respectively, with the to, -1 maximums and minimums of v™. The integer to, which 
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Figure 2. Population fractions Xi for n = 200 and a = 0.15. In black results 
from a simulation after t generations and in gray exp [ -Amv'" t]. (A) and (B): 
Species evenly spaced along the niche axis for t = 1000 and t = 10,000 generations, 
respectively. (C) and (D): Species randomly distributed along the niche axis for 
t = 1000 and t = 10,000. 



II The gray lines, generated from v'", are actually lines. They were drawn tick just to show their 
coincidence with the black thin lines produced by simulations. 
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gives the minimal eigenvalue, is a function of the width a of the niche, m = m{a). It 
does not depend from n provided n is large enough. Nevertheless, as we will show in 
section IV, m becomes a function of n and a for small values of both these parameters. 
For example, for (T=0.15, m-1 = 4 for all even n greater or equal than 8. This lower n 
limit arises because the maximum possible number of peaks that can be accommodated 
with n vector components is n/2 (one half of the components of v™ pointing up and 
the other half down). So in this particular case n/2 must be greater or equal than 4, 
and, in general, n/2 must be greater or equal than m-1. 

Another remarkable result about m is that it is always an odd number (and then 
the number of clumps is even). The reason for this can be traced from the cosines 
appearing in (|T2|) making contributions to the eigenvalues of opposite signs: positive 
for odd k and negative for even k. As a consequence the number of peaks, equal to 
m-1, is always even. 



3. SELF-ORGANIZED SIMILARITY PERSISTS UNDER MORE 
REALISTIC ASSUMPTIONS 



In order to consider more realistic assumptions, abandoning the simplifying conditions 
S1-S4, we rely in the following exclusively on simulations. Since the emphasis in SOS 
is on transient maintenance of clumps of similar species, one might wonder about how 
initial conditions determine the results, and how species that are being driven extinct 
ever managed to get up to high density in the first place. So, as before, the ODE 
system is integrated starting from initial Xi which are random numbers between O and 
1. We checked in all the cases that changes in the initial populations don't introduce 
qualitative changes. 

1. From evenly to randomly distributed species. 

What happens in the general case of randomly distributed species over the niche 
axis? In this case the spectrum and v"* are obtained numerically from A. It turns 
out that simulations produce quite the same results. We illustrate this in Figl2] where 

-, for the particular 



we plot the population fractions normalized to one, Xi = ^^rn 

parameter values n—200 and (t=0.15. The resemblance is clear when comparing panels 
(C) and (D) with, respectively, (A) and (B). In fact, the spectrum of eigenvalues in 
both cases is very similar as it is shown in Table 1 (the values on the right correspond 
to averages among simulations). 





Evenly spaced 


Randomly distributed 


Ai 

A2 

Al98t/i 
Al99t/i 
^200th 


-0.3938 
-0.3938 

45.391 
45.391 
104.387 


-0.4 ±0.01 
-0.4 ±0.01 

46±0.96 

46±0.96 

105 ±1.93 



Table 1: eigenvalues of A for n = 200 & a 

i.e. Am = Ai = 



= 0.15 ordered from small to large 

A2. 



2. Taking into account border effects in a linear niche axis. 

We also analyzed what happens when a linear, instead of a circular niche axis 
(PBC), of length L is considered. The competition coefficients for these open boundary 
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conditions (OBC) are now given by 

. Mi-Mj -2 erf( 

/J ■ ■ ^ p V 2tT f 



2L— /^^— /^j ~ 
2ct - 



crf( 



Mi+Mj 



erf( 



L—fL, 



crf(^; 



(13) 



When using competition coefficients given by ([T^ witli L=l, again, a lumpy pattern 
emerges although it shows some quantitative differences. For example, a four lump 
pattern occurs for smaller values of a, e.g. (j= 0.12 instead of a— 0.15 (panel (A) in 
Figin]). Additionally, although A™ is still negative, due to the factor multiplying the 
Gaussian in (|13p . the matrix A is no longer symmetric and so there appear complex 
eigenvalues. 
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Figure 3. Fraction of species Xi (black bars, left vertical axis) for n = 200, cr 
= 0.1 and open boundary conditions (coefficients given by II13I I writh L=l) after 
500 generations and the corresponding entropy for each lump and gap region Sr 
(gray dashed lines, right vertical axis). (A): Uniform maximum growth rate r and 
carrying capacity K. (B): Varying r (the r^ are random numbers with average 
value equal to 1) and uniform K. (C): Varying r (the r; are random numbers 
with average value equal to 1) and K ( SKmax/K = 0.2, see text) from species 
to species. 



3. The effect of a non uniform growth rate. 
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Simplification S3 was to consider a uniform r. Indeed it is simple to realize that 
an r varying from species to species does not introduce major changes. This is because 
what is relevant for the equilibrium values x* are the terms between brackets in the 
LVCM equations (see panel (B) in Figl3]). 

4-. The effect of the heterogeneity in the carrying capacity. 

We find that when variations ±6Ki of the carrying capacity around an average 
value K occur in such a way that the amplitude of these fluctuations, SK^ax, is 
no greater than 10 % of -R' the lumpy pattern changes but is similar to the one 
corresponding to the homogeneous case. If, in addition, one assumes that the carrying 
capacity of neighbor species along the niche axis have similar carrying capacities and 
larger variations are only possible for species which are far away on the niche axis, 
then larger values of SK^ax/K still preserve SOS (panel (C) in FiglS]). On the other 
hand, strong random variations of the carrying capacity along the niche axis in general 
destroy the SOS pattern. 

In Fig|3] we show the population fractions obtained when the more realistic 
conditions 2 to 4 are gradually taken into account. In the three panels we plot 
the results produced by simulations starting from the same initial distribution of 
populations. Panel (A) corresponds to OBC and homogeneous r and K, panel (B) to 
OBC, heterogeneous r and homogeneous K and panel (C) to OBC and heterogeneous 
r as well as K. Notice that although the lumpy structure becomes less clear as the 
original restrictions are lifted, it is still recognizable in panel (C). In order to provide 
a more quantitative test for the clumping to the favorable niches, it is necessary to 
introduce an observable which measures species coexistence or diversity. Among the 
different indices proposed to measure species diversity perhaps the most common is 
the Shannon- Wiener index [13 ] .114 ) . or in the physics language the well known entropy 
iS*, defined by 

n 

S = — 2_. ^i 111 ^i ■ 
i=\ 

Moreover, entropy analysis has been used to quantify species diversity and niche 
breadth Xh and to recognize ecological structures (see dH and references therein). 
Therefore we proceed as follows. From the homogeneous r and K situation we obtain 
the modulation along the niche axis determining the number and positions of lumps 
and gaps. In this specific case there are four lumps separated by gaps all of the same 
length. Thus we divide the niche axis into 9 regions: 4 lumps and 3 gaps between them, 
all the 7 of length 0.125, plus the two smaller gaps at the niche borders completing 
the remaining length of 0.125. The amount of entropy S^ calculated for each region 
(r=l,2,...), measures the species diversity (represented by gray dashed lines in panels 
(A)-(C)). Notice that the profile of Sr for the three situations is similar although, as 
expected, it offers more clear cut evidence of lumps and gaps for the homogeneous 
situation of panel (A): the entropy is in general lower (higher) at the gaps (lumps of 
coexistence) than in panels (B) or (C). Therefore, we conclude that the considered 
simplifications don't introduce substantial changes and that SOS survives in more 
realistic conditions. 

Other competition kernels and multidimensional niches 

It was argued that the formula ([2]) is a special case and that competition 
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coefficients are typically non-Gaussian [SJ [T^. Some recent analyses explore more 
general non-Gaussian competition kernels of the form |18) : 



= e-(^ 



'.)P 



(14) 



which reduces to the Gaussian one for p=2. Moreover, it was claimed that Gaussian 
competition does not lead to patterns but is a borderline case between patterns and 
non-patterns regimes [lOj . However, this depends on whether or not one takes into 
account the finiteness of the niche axis. When it is taken into account, as we do by 
using a truncated kernel, p—2 is no longer a border case. Rather the lumpy pattern 
occurs for any real kernel exponent p above 1, for example, p — 1.5 as is it is shown 
in panel (A) of Fig. |4]for a— 0.19. This is because the only change in the formula 

for the eigenvalues (|10l) is in the coefficients Cj{n,a) which, for a general value of 

-(L-L)p 



the exponent p, are given by Cj{n,a) 



while the expression pT|) for the 



eigenvectors remains unchanged. Panel (B) of this figure is a plot of the components 
of the V™ showing that its peaks (valleys) coincide with the lumps (gaps). Another 



0.05 




Figure 4. Results for a non-Gaussian icernel with p=1.5, n = 200, cr = 0.19 and 
PBC. (A): Distribution of species for after 2500 generations. (B): Components of 
the dominant eigenvector v™. 



common criticism is that it is not very realistic to consider a one-dimensional niche, 
rather niches (utilizations) in general are multi-dimensional (THl [ID] It turns out that a 
multi-dimensional niche only makes the math a little bit less straightforward. Suppose 
that the n species are distributed at random in a 2-dimensional niche with axes fii 
and ^2- Then one can assign an index i to each population, located at the point in 
niche space given by a couple (/.iii,/i2i)i ^^'^ group them into a vector of n components. 
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Therefore, the expression for the competitfon coefficient between species i, focated in 
this niche space at a point of coordinates {pii,fi2i), and species j, at {fiij^fj,2j), can be 
written as 

aij=e (^^ . (15) 

It turns out that this preserves the cyclic property of the A matrix - its rows are 
cychc permutations of the first one -, a property required to get the expressions for 
the eigenvalues (fTO|) and the eigenvectors (fTTj) [12]. Fig. [5]shows the results for a— 0.2. 
It shows a general result we found: if in the case of a one-dimensional niche {d—1), 
for a given value of a, there are m-1 lumps, for a two-dimensional niche id=2) there 
occur (m-1) X (m-1) lumps (for a = 0.2 there are 2 lumps for d—1 while for d—2 there 
are 4 lumps). 




Figure 5. Distribution of species in a two-dimensional niciie of coordinate axis 
Hi and 112 for n = 15x15= 225 species, and a = 0.2, PBC after 100 generations. 
2x2=4 clumps are observable. 
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4. THE DEPENDENCE OF CLUMPING ON THE NICHE WIDTH 
AND THE CLUMPING TRANSITION 

How close species can be packed along the niche axis is commonly measured by the 
parameter d/a, where d is the separation between species jl . In the case of the model 
under consideration, d can either measure the separation between a) lumped groups 
of species, persisting during long transients or b) surviving species (one per lump), 
for asymptotic times. So we get either an estimate for the species packing or for the 
group packing. In any event, this distance coincides with the inverse of the number of 
peaks of v™, which depends on a, and is given by n^o (cr) = m{cr)-l. FiglHlshows rioo 
for a ranging from 0.05 to 0.5 and the number of species fixed to n = 200. There is a 
series of steps, located at values a-g, that become wider as as increases. The height of 
these steps is always 2. This is because, as we have seen in section 2, the number of 
peaks of v*" is always an even number. That is, if ctJ (tr^) corresponds to a tending 
to (Ts by the left (right), then noo(o'7) = noo(CT+)+2. For example, if cr > Cs ~ 0.169 
then V™ has always 2 peaks, below this value the number of its peaks jumps to 4, 
and so on. We find that when a tends to af , rioo can be fitted with the power-law 
0.09a~^'^^ (dashed line in Fig|6]). The packing parameter in the different step regions 
can be approximated, by taking the number of peaks at each as as the semi-sum of 
the numbers of peaks at each side of the step, as: 



d 1 1 



+ ^_L^ l^-\\ ^ I nnn^-0.75- 



as asl/2{noo{al) + noc{as )) cts + 0.09cr, 



(16) 



This quotient varies from approximately 1.96 for the first step, at as ~0.169, to 1.1 for 
the last step, at as — 0.05. This is in good agreement with many field measurements 
that found a species packing ratio always lying between 1 and 2 [1]. It is worth 
remarking that when a decreases, Am -which is in general negative- increases, until at 
some critical value, CTc, it becomes 0. That is, in Thom's catastrophe theory language 
|21| . a degenerate critical point or a Non-Morse critical point. This Cc depends on 
the number of species: it decreases with n. We computed adn) as the values such 
that Xm{n,ac) becomes 0. In Figl7]we show this. Notice that for n > 40 Uc scales 
as 2n~^, i.e. the double of the initial average separation between species. As a is 
decreased in simulations -for a fixed value of n- so that it becomes closer and closer 
to ac and Am moves towards 0, we observed that the time to reach the lumpy pattern 
grows unbounded. This is the well known phenomenon of critical slowing down^^ : 
the characteristic relaxation time of the dominant eigenmode is proportional to 1/A„i- 
In fact, taking ?i=200, for a = 0.15 the 4 lumps are noticeable after typically 500 
generations while the 6 lumps for a = 0.1 require around 20,000 generations and for 
a = 0.075 it takes a huge number of generations (more than 500,000) to produce the 
10 lumps pattern. The inset of FigUis a zoom of A™ vs. a. It shows that, at least 
for all practical purposes, the clumping becomes noticeable at cr ~ 0.075. 

5. CONCLUSIONS AND FINAL COMMENTS 

A realistic mathematical description of the dynamics of a large number of species 
placed along a resource spectrum is a complicated issue for which an exact solution 
is not available. In fact, analytical work looks at the long-term equilibria of models. 
The alternative to deal with the transients are simulations. 
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Figure 6. Number peaks rioo of the dominant eigenvector v™ as a function of cr 
(ra = 200). The jumps follow a power-law distribution indicated by a dashed line. 



However a simulation approach, like the one used by Scheffer and van Nes [3] 
may leave room for doubts on whether things might be artifacts. We made a series of 
simplifications which allow an analytic proof, by working directly with the community 
matrix A, of the emergence of SOS. Roughly, the lumpy pattern one is seeing is the 
exponential of the dominant eigenvector v™ of A (multiplied by the time). 

We later have shown that this is indeed a robust result. The clumping 
phenomenon does not depend on the boundary conditions nor on the kernel exponent 
p (provided it is greater than one), it is quite independent on the heterogeneity of 
species parameters r and K, it occurs for a wide range of a, and in more than one 
niche dimensions. In addition, Roelke and Eldridge [23j found similar patterns in a 
different resource competition model. They also suggest that this mechanism is not 
very fragile. Additional supporting evidences can be found in p4l. 

A crucial element to get clumping is to take into account the finiteness of the niche 
axis. This, besides being realistic, leads to a A™ with a negative real part (and then to 
lumps and gaps). We want to remark that either OBC or the standard implementation 
of PBC imply a niche which is finite. This explains remarkable differences with the 
outcomes reported in [10]. The procedure they use to implement PBC consists in 
taking a periodic array of copies of the same system. This "perfectly periodic" 
boundary conditions mimic an infinite niche axis. The first of such differences is 
that, in our case, the SOS is robust against variations on the kernel: a negative Am is 
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Figure 7. Log-log plot of ac vs. n for n between 20 and 200. The dashed line 
corresponds to 2/n. 



obtained whenever the exponent p of the kernel interaction is a real number greater 
than 1. The second difference, is that the parameter a, controlhng the width of 
each species distribution, plays a fundamental role. There is a critical value, ac, below 
which there is no clustering. On the other hand, in a virtually infinite niche axis, since 
it is always possible to set ct = 1 by rescaling fj,, the clustering should not depend on 
a. However it seems natural that things in ecosystems depend on a, and usually the 
interest is precisely in measuring this effect. This is another powerful reason to prefer 
an implementation of boundary conditions like the one we are using. 

Similar approaches to determine pattern formation in phenotype space have also 
been used by Levin and Segel [25], Sasaki p6l and more recently by Meszena and 
co-workers [22]j[2H]- The main difference of our approach is that a discrete set of 
phenotypes is considered, instead of a continuum. This is an important feature, since 
firstly it can affect the assessment concerning how robust SOS is. That is, in the case of 
a continuous set of phenotypes, it was shown that an arbitrarily small perturbation can 
destroy the continuous coexistence making the species distribution discrete [27].[28]. 
And this could be taken as an evidence of the breaking of SOS. On the other hand, 
in our model, from the very beginning, the distribution of species is discrete and SOS 
it is rather understood as a strong overlap of the species distributions. Secondly, it 
allows to go beyond the large n limit and applying it to real communities, involving a 
number of species n of intermediate size. For example, this was tested for the case of 
phytoplankton communities in a lake ecosystem involving between 50 to 100 species 
and the agreement between theory and empirical data is quite good [29]. Moreover, 
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in this real ecosystem the fact that a is not the same for all species, rather varies from 
species to species, doesn't spoil the lumpy pattern. 

There are interesting parallels with similar phenomena in physical systems: 

• For example, the fact that the emergence of the lumpy pattern is related to the 
eigenvector with the minimum eigenvalue and the number of lumps is determined 
mostly by the model parameter a resembles the spinodal decomposition|3l That is, 
under the spinodal decomposition the system develops a spatially modulated order 
parameter whose amplitude grows continuously from zero and extend throughout 
the entire system. This results in domains of a characteristic length scale called 
the spinodal length Asp which usually depends strongly on temperature (because 
the second derivative of the free energy becomes increasingly negative deep inside 
the region delimited by the spinodal) [50] . 

• The emergence of power laws and critical slowing down are attractive ingredients 
to physicists since they are signatures of self-organized criticality (SOC) [31] . 
This tendency to spontaneously self-organize into a critical state, without any 
significant "tuning" of some control parameter, usually refiects a share of the 
same fundamental dynamics for many different systems referred to as universality. 
While the origin of critical slowing down is clear explained by the existence of a 
degenerate critical point, the power law distribution for the plateaus of the number 
of lumps vs. a is not completely understood and deserves further analysis. 

• The application of techniques and concepts of statistical mechanics into very 
different realms like ecosystems might be of interest to ecologists, statistical 
physicists and to the growing community on the intersection of both fields. In 
that sense, the analytical proof of clumping is based on statistical mechanics 
results from Berlin and Kac [12] when they were analysing the spherical model 
of a ferromagnet. 

• The calculation of the entropy for different regions of a system, to get an overview 
about the level of correlation between elements in each region has been used 
in several contexts closer to physics. For example: cellular automata |32) . 
deterministic models of nonlinear dynamics j33j, glass-forming materials [34], 
astrophysics of galaxies and clusters [35], image processing [36], to mention some. 
In our case it has shown to be useful to identify the lumpy structure. 

To conclude, it is remarkable that the predictions on the number of groups of 
species that can be packed along the niche axis are quantitatively consistent with field 
data for a wide range of values of both the width of the niche and the number of 
species. 
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